#####################
# DID panel data
#####################

rm (list = ls())

setwd("~/Dropbox/Climate change and disasters/07_replication")

#sink("04_analysis_climate_covariates_partisanship.txt")

library(foreign)
library(stargazer)
library(Hmisc)
library(ggplot2)
library(tidyverse)
library(readstata13)
library(did)

###############################
# Load and prepare clean data
###############################

# load
load("disasters_attitudes.RData")
dim(d)

# partisanship
table(d$pid3_10, exclude =NULL)
d$democrat = 0
d$democrat[d$pid3_10==1] =1
d$republican = 0
d$republican[d$pid3_10==2] =1
table(d$democrat, exclude = NULL)
table(d$republican, exclude = NULL)

#################################################
# Outcome: descending; Treatment: previous year
################################################

# fire
set.seed(11)
out1 <- att_gt(yname = "outcome_descending",
              gname = "first_treat_fire",
              idname = "id",
              tname = "survey",
              xformla = ~democrat + republican,
              data = d)
es1 <- aggte(out1, type = "dynamic")
summary(es1)

pe1 <- es1$overall.att
se1 <- es1$overall.se
results1 <- rbind(pe1,se1) # add baselines
results1

event_study_diff_never1 <-   data.frame(
  type          = "dynamic",
  term = paste0('ATT(', es1$egt, ")"),
  event.time= es1$egt,
  estimate  = es1$att.egt,
  std.error = es1$se.egt,
  conf.low  = es1$att.egt - es1$crit.val.egt * es1$se.egt,
  conf.high = es1$att.egt + es1$crit.val.egt  * es1$se.egt)

p_es1 <- ggplot(data = event_study_diff_never1,
                      mapping = aes(x = event.time, y = estimate)) +
  geom_line(size = 0.5, alpha = 2, colour = "black") +
  theme_light()  +
  geom_hline(yintercept = 0, colour="black",  linetype = "dotted")+
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), show.legend = FALSE, linetype= 1, size = 1.1,
                  color=ifelse(event_study_diff_never1$event.time<0,"gray60","black"))+
  #geom_pointrange(aes(ymin = point.conf.low, ymax = point.conf.high), show.legend = FALSE, size = 1.1)+
  xlab("Panel wave since exposure") +
  ylab("DiD Estimate") +
  scale_x_continuous(breaks = c(-1:1)) +
  scale_y_continuous(limits = c(-0.3,0.3)) +
  theme(axis.text.y = element_text(size = 14),axis.title.y = element_text(size = 16)) +
  theme(axis.text.x = element_text(size = 14),axis.title.x = element_text(size = 16)) +
  theme(axis.title = element_text(color="black",  size = 12))+
  annotate(geom="text", x=0.85, y=7.6, label="",
           color="black",
           size = 8)
f1 = p_es1 +  
  ggtitle("Effect of exposure to fires") + 
  theme(plot.title = element_text(size = 18)) +
  theme(plot.title = element_text(hjust = 0.5)) 

# flood
set.seed(11)
out2 <- att_gt(yname = "outcome_descending",
              gname = "first_treat_flood",
              idname = "id",
              tname = "survey",
              xformla = ~democrat + republican,
              data = d)
es2 <- aggte(out2, type = "dynamic")
summary(es2)

pe2 <- es2$overall.att
se2 <- es2$overall.se
results2 <- rbind(pe2,se2) # add baselines
results2

event_study_diff_never2 <-   data.frame(
  type          = "dynamic",
  term = paste0('ATT(', es2$egt, ")"),
  event.time= es2$egt,
  estimate  = es2$att.egt,
  std.error = es2$se.egt,
  conf.low  = es2$att.egt - es2$crit.val.egt * es2$se.egt,
  conf.high = es2$att.egt + es2$crit.val.egt  * es2$se.egt)

p_es2 <- ggplot(data = event_study_diff_never2,
                mapping = aes(x = event.time, y = estimate)) +
  geom_line(size = 0.5, alpha = 2, colour = "black") +
  theme_light()  +
  geom_hline(yintercept = 0, colour="black",  linetype = "dotted")+
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), show.legend = FALSE, linetype= 1, size = 1.1,
                  color=ifelse(event_study_diff_never2$event.time<0,"gray60","black"))+
  #geom_pointrange(aes(ymin = point.conf.low, ymax = point.conf.high), show.legend = FALSE, size = 1.1)+
  xlab("Panel wave since exposure") +
  ylab("DiD Estimate") +
  scale_x_continuous(breaks = c(-1:1)) +
  scale_y_continuous(limits = c(-0.3,0.3)) +
  theme(axis.text.y = element_text(size = 14),axis.title.y = element_text(size = 16)) +
  theme(axis.text.x = element_text(size = 14),axis.title.x = element_text(size = 16)) +
  annotate(geom="text", x=0.85, y=7.6, label="",
           color="black",
           size = 8)
f2 = p_es2 +  
  ggtitle("Effect of exposure to floods") + 
  theme(plot.title = element_text(size = 18)) +
  theme(plot.title = element_text(hjust = 0.5)) 

# severestorm
set.seed(11)
out3 <- att_gt(yname = "outcome_descending",
              gname = "first_treat_severestorm",
              idname = "id",
              tname = "survey",
              xformla = ~democrat + republican,
              alp = 0.1,
              #control_group = c("notyettreated"),
              data = d)
es3 <- aggte(out3, type = "dynamic")
summary(es3)

pe3 <- es3$overall.att
se3 <- es3$overall.se
results3 <- rbind(pe3,se3) # add baselines
results3

event_study_diff_never3 <-   data.frame(
  type          = "dynamic",
  term = paste0('ATT(', es3$egt, ")"),
  event.time= es3$egt,
  estimate  = es3$att.egt,
  std.error = es3$se.egt,
  conf.low  = es3$att.egt - es3$crit.val.egt * es3$se.egt,
  conf.high = es3$att.egt + es3$crit.val.egt  * es3$se.egt)

p_es3 <- ggplot(data = event_study_diff_never3,
                mapping = aes(x = event.time, y = estimate)) +
  geom_line(size = 0.5, alpha = 2, colour = "black") +
  theme_light()  +
  geom_hline(yintercept = 0, colour="black",  linetype = "dotted")+
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), show.legend = FALSE, linetype= 1, size = 1.1,
                  color=ifelse(event_study_diff_never3$event.time<0,"gray60","black"))+
  #geom_pointrange(aes(ymin = point.conf.low, ymax = point.conf.high), show.legend = FALSE, size = 1.1)+
  xlab("Panel wave since exposure") +
  ylab("DiD Estimate") +
  scale_x_continuous(breaks = c(-1:1)) +
  scale_y_continuous(limits = c(-0.3,0.3)) +
  theme(axis.text.y = element_text(size = 14),axis.title.y = element_text(size = 16)) +
  theme(axis.text.x = element_text(size = 14),axis.title.x = element_text(size = 16)) +
  annotate(geom="text", x=0.85, y=7.6, label="",
           color="black",
           size = 8)

f3 = p_es3 +  
  ggtitle("Effect of exposure to severe storms") + 
  theme(plot.title = element_text(size = 18)) +
  theme(plot.title = element_text(hjust = 0.5))

# hurricane
set.seed(11)
out4 <- att_gt(yname = "outcome_descending",
              gname = "first_treat_hurricane",
              idname = "id",
              tname = "survey",
              xformla = ~democrat + republican,
              alp = 0.1,
              #control_group = c("notyettreated"),
              data = d)
es4 <- aggte(out4, type = "dynamic")
summary(es4)

event_study_diff_never4 <-   data.frame(
  type          = "dynamic",
  term = paste0('ATT(', es4$egt, ")"),
  event.time= es4$egt,
  estimate  = es4$att.egt,
  std.error = es4$se.egt,
  conf.low  = es4$att.egt - es4$crit.val.egt * es4$se.egt,
  conf.high = es4$att.egt + es4$crit.val.egt  * es4$se.egt)

p_es4 <- ggplot(data = event_study_diff_never4,
                mapping = aes(x = event.time, y = estimate)) +
  geom_line(size = 0.5, alpha = 2, colour = "black") +
  theme_light()  +
  geom_hline(yintercept = 0, colour="black",  linetype = "dotted")+
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), show.legend = FALSE, linetype= 1, size = 1.1,
                  color=ifelse(event_study_diff_never4$event.time<0,"gray60","black"))+
  #geom_pointrange(aes(ymin = point.conf.low, ymax = point.conf.high), show.legend = FALSE, size = 1.1)+
  xlab("Panel wave since exposure") +
  ylab("DiD Estimate") +
  scale_x_continuous(breaks = c(-1:1)) +
  scale_y_continuous(limits = c(-0.3,0.3)) +
  theme(axis.text.y = element_text(size = 14),axis.title.y = element_text(size = 16)) +
  theme(axis.text.x = element_text(size = 14),axis.title.x = element_text(size = 16)) +
  annotate(geom="text", x=0.85, y=7.6, label="",
           color="black",
           size = 8)
f4 = p_es4 +  
  ggtitle("Effect of exposure to hurricanes") + 
  theme(plot.title = element_text(size = 18)) +
  theme(plot.title = element_text(hjust = 0.5)) 

###############
# Final plot
###############

library(ggplot2)
library(gridExtra)

tiff('FigA2.tiff', units="in", width=10, height=7, res=200, compression = 'lzw')
grid.arrange(f1,f2,f3,f4, ncol = 2, top = "")
dev.off()

sink()

